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The physical properties of granular materials have been extensively studied in recent years. So 
far, however, there exists no theoretical framework which can explain the observations in a unified 
manner beyond the phenomenological jamming diagram This work focuses on the case of static 
granular matter, where we have constructed a statistical ensemble 01 which mirrors equilibrium sta- 
tistical mechanics. This ensemble, which is based on the conservation properties of the stress tensor, 
is distinct from the original Edwards ensemble and applies to packings of deformable grains. We 
combine it with a field theoretical analysis of the packings, where the field is the Airy stress function 
derived from the force and torque balance conditions. In this framework. Point J characterized by 
a diverging stiffness of the pressure fluctuations. Separately, we present a phenomenological mean- 
fleld theory of the jamming transition, which incorporates the mean contact number as a variable. 
We link both approaches in the context of the marginal rigidity picture proposed by 0, U • 



I. INTRODUCTION 



B. Towards a non-equilibrium theory 



A. Granular media 

Granular materials are ubiquitious in everyday life, as 
sand, salt or rice; or as the breakfast cereal in your bowl. 
They are also of fundamental importance in industrial 
processing, from the transport of coal to drug manufac- 
turing. The shape of the surface of the Earth is largely 
determined by the erosion of rock to sand and dust and 
the patterns they form subsequently. If we want to truly 
understand the world we live in, understanding granular 
materials must be part of our efforts. Yet, our grasp of 
the behavior of granular matter is limited [a, S [I]- 
A consensus about what the important questions in the 
field arc, and on suitable approaches to answer them, has 
only started to emerge in the last few years [l|, . 

The chief obstacle that hinders our understanding is 
that granular matter is fundamentally out of thermal 
equilibrium, since the typical interaction energies of the 
particles are many orders of magnitude larger than k{,T. 
Therefore, granular materials are effectively at T = 
from a conventional point of view. This means that in 
equilibrium, the system does not explore the space of 
states available to it, thus leading to a breakdown of er- 
godicity. It also implies that the Boltzmann ensemble is 
irrelevant, and we are left without the framework of equi- 
librium statistical mechanics to guide our understanding. 

Since granular systems do not equilibrate sponta- 
ncusly, the method of preparation of a granular packing 
matters, i.e. granular materials are history-dependent. 
A steady-state can only be achieved through highly con- 
trolled experimental conditions like repeated tapping [Tol| 
or compression This complicates the study of gran- 
ular materials since it is far from clear if a universal fram- 
work is even possible. 
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Despite the obstacles described in the previous section, 
some progress has been made towards a non-equilibrium 
theory of granular materials. Most theoretical work has 
attempted to extend concepts like temperature or en- 
tropy from the thermodynamic context to the granular 
context. 

It has been recognized that granular packings and 
granular systems under shear or other forms of agita- 
tion share several important properties of glassy systems. 
Like granular materials, glassy materials have fallen out 
of equilibrium in the sense that thermal fluctuations are 
too small to allow the system to undergo rearrangements 
on short time scales. Some model glass formers, like 
Lennard-Jones fluids, are structurally similar to granu- 
lar packings. A common concept that has emerged for 
both systems is the idea of a rugged potential energy 
landscape with a myriad of energy minima. For glassy 
systems, these are termed the inherent structures 12| in 
which the system spends a long time before escaping over 
an energy barrier to a different one via an activated pro- 
cess. For granular systems, the potential energy minima 
correspond to packings in mechanical equilibrium, and 
the dynamics are provided by tapping or slow shearing. 
For situations where the deformation is gradual enough 
for the system to move by going through a sequence 
of separate rearrangements (the quasistatic shear limit), 
Bouchaud's trap model for glasses [13] has been adapted 
to include shearing as the SGR (Soft Glassy Rheology) 
model by Sollich [ij. 

The liquid to solid phase transition in granular media is 
also known as the jamming transition. It is controlled by 
parameters including an external driving force like shak- 
ing or shearing, but that do not include temperature. 
This transition shares some of the properties of the glass 
transition. In the glass transition the system changes 
from a liquid to a disordered solid without an obvious 
change in microscropic structure, unlike for the conven- 
tional liquid-solid transition where the system goes from 
a state of high symmetry (the liquid phase) to a state of 
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lower, broken symmetry (the crystalline solid) p^ . 

These similarities have prompted Liu and Nagel to pro- 
pose [i| a unified phase diagram for glassy and granular 
materials, as well as related materials like foams, bubbles 
and colloids. Glassy materials are in the temperature - 
density plane of the diagram. The transition between 
a granular solid and a granular liquid with increasing 
shear is in the density-shear plane. The packing fi'action 
axis shared by both types of materials is characterized 
by granular packings in mechanical equilibrium, or alter- 
natively glassy systems at zero temperature. 

The properties of the phase transition depend on how 
the phase boundary is approached. Along the packing 
fraction axis, reducing the density of an isotropic jammed 
assembly of grains leads to a point at which both the bulk 
modulus and the shear modulus, i.e. the resistance of the 
system to uniform compression and shear, vanish 
The transition point has been termed Point J by [l| , and 
it has some of the features of a critical point, though it 
not obvious that Point J has any bearing on the glass 
transition [l8l |. 

O'Hern et al. [13] have found power-law scalings and 
apparently universal properties when J is approached 
from the high-density limit, and scaling as well as a diver- 
gent length scale have also been observed by Drocco et al. 
[3] and Teitel et al. [13] when J is approached from the 
low-density limit. Experimentally, the microscopic prop- 
erties of static packings, such as the force and contact 
number distributions, as Point J is aprroached along the 
packing fraction axis have been studied by Majmudar et 
al. lU. 

For spherical particles, Point J can be linked to the 
underlying microscopic properties of the granular pack- 
ings. For spherical granular materials with purely repul- 
sive, short-range interactions. Point J is identical to the 
isostatic point, the point where the number of degrees of 
freedom is equal to the number of constraints imposed by 
the intergranular forces 17]. In the frictionless case, the 
jamming thresholds of different configurations converge 
with increasing system size to a single packing fraction 
[l7| which is identical to random close packing (RCP), 
the highest packing fraction that can be attained by ran- 
dom packings of hard grains [2^ . For three dimensional 
packings of spheres, 4>j = 0.63, while for two dimensional 
packings of disks it is given by 0,/ = 0.84. 

A majority of experimental studies of the jamming 
transition probe the transition as a function of the shear 
rate, where the phase boundary is crossed at finite posi- 
tive pressure. In this case, only the shear modulus van- 
ishes at the transition and the nature of the transition is 
different from the compression case. The signature of this 
transition is intermittent jamming and unjamming with 
large-scale stress fluctuations, and pronounced dilatancy 
effects [§]. A different type of jamming occurs in systems 
where there is a sustained motion of the grains, and ki- 
netic energy plays a role. These transitions are charac- 
terized by the appearence of large scale spatio-temporal 
fluctuations in the particle motion j23l . [24 ] . These fluc- 



tuations have been termed dynamical heterogeneities, af- 
ter the very similar behavior which is observed in super- 
cooled liquids [1^ . Further analysis using tools developed 
for the glass transition, like the self-intermediate scatter- 
ing function, has revealed a growing l eng th scale as the 
jamming transition is approached [2a. |27| |28|. 

C. Granular statistical mechanics 

Going further in the analogy to thermodynamics and 
ordinary statistical mechanics, several attempts have 
been made to build a framework equivalent to equilib- 
rium statistical mechanics for granular materials. 

The first attempt was by Edwards [29!] who replaced 
the energy conservation law for thermal systems with a 
volume conservation law for granular packings. Edwards 
makes the microcanonical assumption that all states at 
the same free volume V are accessed with the same proba- 
bility. For a subsystem with volume v of the full packing, 
one then obtains a canonical distribution for the proba- 
bility to access the state 

= exp(-t;,/x(V^)); (1) 

where the Lagrange multiplier x, or the compactivity, 
plays the role of temperature, x is the best known ex- 
ample of a granular temperature for static packings (see 
[20] for a review). 

The validity of the Edwards ensemble has been inves- 
tigated for experimental and simulated grain packings. 
Some studies have been direct investigations of the free 
volume distribution in packings [U, [H, [H, [sj. A lot 
of the tests have been indirect, generally through the 
fluctuation-dissipation theorem [1^ (see [3a] and [Sq). 

Other out-of-equilibrium temperatures , mostly based 
on the fluctuation-dissipation theorem, have been pro- 
posed in the context of dynamical granular media. One 
study 37] measured several of them in a simulated sys- 
tem of a two-dimensional sheared foam, and found that 
the different temperatures agree with each other, show- 
ing that effective temperatures are a useful concept for 
driven granular systems. 

In recent years, it has become more common to apply 
the Edwards equiprobability hypothesis to force config- 
urations on static granular packings (see e.g. [11], [s^] 
and [40] ) . An equivalent approach to the case of the vol- 
ume can be taken, and a force-based canonical ensemble 
is then introduced. The resulting granular temperature 
has been dubbed the angoricity by Edwards [4l| . 

D. This work 

Here we derive the conservation law responsible for the 
force ensemble from first principles, and show the condi- 
tions under which the canonical ensemble is valid (see 
also our previous work [2]). We construct and exam- 
ine a statistical mechanics-like framework in which Point 
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J and the properties of jammed systems can be studied. 
We restrict our attention to packings on the density-shear 
plane of the jamming phase diagram, where we find a con- 
servation law upon which we base an ensemble in which 
stress plays the role of energy. We then use the basis 
laid out by the ensemble to develop a field theoretical 
model for jammed packings, which predicts a diverging 
stiffness at the jamming transition. Finally, we discuss 
an empirical mean-field theory derived from the analysis 
of simulated packings. 

Section 2 derives a conservation law for the force- 
moment tensor, i.e. the volume integral of the stress ten- 
sor, from the boundary property shown in [3, \^^. This 
forms the basis for the derivation of the stress ensemble, 
where the force moment tensor plays the role of energy, 
from first principles. We derive a tensorial analog of tem- 
perature, the angoricity, which is related to stress fluc- 
tuations rather than energy. The concept of angoricity 
has been introduced earlier. To avoid any confusion, we 
precisely define our usage of the term in this paper. We 
also discuss the case of an isotropically compressed sys- 
tem where we can simplify the ensemble to a scalar stress 
variable and a scalar temperature. The stress ensemble 
has been numerically tested in Q , and we provide a sum- 
mary of the results. Finally, we discuss an exact model 
of the jamming transition compatible with the universal 
equation of state found in the simulation. 

Section 3 builds a field theoretical model of granular 
packings based on the scalar stress function appropriate 
for isotropic compression. Based on symmetry arguments 
and on comparison to simulation results, we find that the 
packings can be described by a laplacian field theory with 
a stiffness that depends on the compression. The stiffness 
diverges at the jamming transition, so that the entropy 
defined by the number of packings available vanishes. Fi- 
nally, we discuss shear fiuctuations in the context of the 
field theory. 

Section 4 combines the equation of state found empir- 
ically with the form of the density of states. It incorpo- 
rates another empirical relation linking stress and contact 
numbers, which has also been derived from a stability ar- 
gument, to define an effective free energy. We study the 
predictions of the effective free energy for the jamming 
transistion, and find that the theory admits divergent 
fluctuations at Point J; but only if formulated in terms 
of a variable u measuring the deviations from the sta- 
bility line. We conclude by showing that the simulated 
system indeed shows a divergence in the fluctuations of 
u as Point J is approached. 



II. THE STRESS ENSEMBLE 

In this section, we derive from flrst principles an en- 
semble based on the conservation of the total force mo- 
ment tensor for static packings of granular materials un- 
der given boundary conditions. This stress ensemble pro- 
vides the formal basis of the force network ensemble used 



in [H, li^l . We have already presented a derivation of the 
ensemble in 0, however the derivation below has been 
improved by taking into account the general formalism 
for nonequilibrium temperatures presented by Bertin et 
al. ,43] , and it has been generalized to the full force mo- 
ment tensor S. 



A. Conserved quantities in stable grain packings 

The total stress tensor of a granular packing with vol- 
ume V and only contact forces can be written as [i^ 

where the sum runs over all the contacts in V, the Fij 
are the contact forces between grains i and j and the 
fij are the distances between grain centers and contact 
points. The extensive quantity S, the total force moment 
tensor, depends only on the boundary conditions of the 
packing which determine a and the total volume V . This 
fundamental result is ultimately a consequence of force 
balance, and it can be demonstrated in a more formal 
way as shown below: 

Ball and Blumenfeld jl^l have shown that local force 
balance in a two-dimensional granular packing can be in- 
corporated through the definition of a vector height field 
field, hfj_, on the loops fi dual to the contact network. 
With this definition, the stress tensor associated with an 
area A can be written as a boundary term that only in- 
volves the boundary vectors c?^j_ and the boundary height 

field hi,, see also Figure 1 in [2]. In three dimensions, the 
same result can be obtained by a continuum approach 
involving a tensorial equivalent of the height field [44] . 

The global force moment tensor is then (the first ex- 
pression is for two-dimensional systems, while the second 
expression is for three dimensional sytems) 

S=2d 4i,^Mb =3d / WxndA. (3) 

bound. "'^^ 

which is the global equivalent to writing the stress tensor 
as the curl of the height field (T = Vx/i((7 = VxM^in 
three dimensions). 

The fact that E can be written as a boundary integral 
implies a conservation principle for force-balanced, static 
grain configurations: each rearrangement of the grains 
within a region of packing, so long as it does not affect 
the boundary of the region, will not change the total force 
moment S of the region. 

The force moment tensor is an extensive quantity, it 
is additive if we neglect boundary effects, as is the case 
for energy in thermodynamics, and it is conserved for any 
dynamics which keeps the boundary conditions of a gran- 
ular packing intact. We now use E to define a canonical 
ensemble, equivalent to the Boltzmann ensemble, where 



4 



it plays the role of the energy (incidentally, E also has 
the dimensions of an energy) . 

An example of rearrangements which conserve the to- 
tal force moment tensor are the wheel moves introduced 
by Tighe et al. on the triangular lattice . 



B. Derivation of the canonical stress ensemble 

Counting the number of states The first step towards 
constructing an ensemble equivalent to the Boltzmann 
ensemble is to define a density of states. The total force 
moment tensor S of a compact area An with TV grains 
is not affected by internal rearrangements. Hence we can 
partition the phase space of all the possible force- and 
torque-balanced packings on An (the blocked states) into 
sectors A/g with the same total force moment tensor E. 
We then count the number of states in each sector AI-^, 
/3Ar(S). This number includes all the possible geomet- 
ric configurations of the grains and will thus always be 
greater than 1, even at the isostatic point. 

In a thermal system, pn{E), the total number of con- 
figurations at a given energy, would form the basis of 
the microcanonical ensemble. A fundamental assump- 
tion underlying thermodynamics and equilibrium statis- 
tical mechanics is that the system is ergodic such that all 
states with a given energy are visited equally. The proba- 
bility of finding a state is then just a flat measure in state 
space, i.e Pv\e — ^/pn{E), and we can then identify the 
density of states as flN{E) = pn[E). 

In a granular system, the dynamical processes which 
move the system from one state to the next are highly 
varied and to a certain degree arbitrary. States within a 
sector can (but may not) be connected by local dynamics; 
but states in different sectors are inaccessible through 
purely local dynamics. We can think of the phase space of 
blocked states through a landscape, similar to the energy 
landscape often invoked for glassy materials. However, 
it is the chosen dynamics, not thermal fluctuations, that 
move the system from one blocked state to another. The 
dynamics are in general non-ergodic, and in the extreme 
case of a purely static system, we recover the T — 0-limit 
of the glass. 

Detailed balance is violated, and hence the system is 
not in thermal equilibrium and the states are not nec- 
essarily accessed with equal probability. Therefore, the 
density of states for Mg, is not only dependent on the 
number of states pAr(I]) at that E, but also on the fre- 
quency Z?!'^^"'' with which each state is accessed by the 
dynamics chosen to create the packing. The quantity 
equivalent to the density of states is then given by: 



(4) 



Preconditions We develop the ensemble along the 
same lines as for a thermal system [45| . Bertin et al. [i^ 
have shown that intensive quantities, which can be inter- 
preted as nonequilibrium temperatures, can be defined 
in a system in steady-state, so long as there is an ad- 
ditive conserved quantity in the system. The conserved 
quantity plays the role of energy, and leads then to a 
temperature variable conjuguate to it. We follow their 
method of derivation below. 

For mechanically stable granular packings, E is the ad- 
ditive quantity that replaces the energy. The following 
results are dynamics-dependent, and stay valid for any 
dynamics which satisfies equation [5l but only as long as 
the same dynamics are used in the preparation of any sys- 
tems we wish to compare (e.g. in steady-state). Whether 
different types of dynamics can give rise to the same den- 
sity of states, and if so, which classes do, is a question 
which we address at the end of this section. 

Consider a system S which contains N grains, with 
total force moment tensor E. We partition the system 
into two compact subsystems 5*1 and S2 with grain 
numbers Ni and and total force moment tensors Ei 
and S2, respectively. Since E is additive, we always have 
S = El -hE2. 

The microcanonical ensemble If we fix the total force 
moment tensor E in S", the conditional probability of find- 
ing a force moment tensor Ei in subsystem Si is defined 
by P(Ei) = P(Ei|E)P(E). Then we can use the def- 
inition of the density of states equation 3] to write the 
conditional probability as 

P(Ei|E) = ^Nity^Y. -Ei)<5(E.,-(E-Ei)). 

If the frequency with which the subsystems are ac- 
cessed factorizes, i.e. if 



p(dyn) _ p(dyn) p{dyn) 

the conditional probability becomes 



P(Si|E) = 



r!Ar,(Ei)»Ar,(E-Ei) 



(5) 



(6) 



as a function of the densities of states of the subsystems. 
Equation [5] is the central assumption in the derivation of 
the stress-ensemble, and it is conceptually equivalent to 
requiring that the dynamics choose state vi of Si inde- 
pendently of state 1^2 of S2 ■ Since the subsystems interact 
through their shared boundary only, we expect the cor- 
rection to equation [5] to scale as 0(l/\//Vi), where A^i 
is the number of grains in 5*1. However, If the system is 
correlated over length scales ^ > VTVi, we expect equa- 
tion [5] to break down as well. We discuss the validity of 
this assumption, and how we have tested it, at the end 
of this section. 

The most probable value E^ at a given E is found by 
setting the derivative of the conditional probability with 
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respect to Si to zero. Since Ei is a tensor, this needs 
to be done for each component separately. We use the 
logarithmic derivative to simplify the calculation, 







01n»jv,(I]i 



Is-- 



91nf7Ar,(i]-Si) 



where we have replaced the derivative in the QN^-teiia 
by d^ij — —d^ij. The first and the second term are then 
opposites of each other, and we define the microcanonical 
equivalent of the inverse temperature by 



9 In r!^, (Si), 51nr!Ar,(E2), 



We will show below that Aij acts indeed like an inverse 
temperature, in that it is independent of the partitioning 
of S*. 



The canonical ensemble To define the canonical en- 
semble, we consider the same system S, but we divide it 
now into one very small partition S'm, with m << and 
the remaining Sn-th (note that we still need m >> 1 for 
the factorization of the density of states to hold). We 
focus our attention now on the small partition where the 
total stress can fluctuate while the full system S with 
fixed E acts as a reservoir, similar to the thermal case. 
We still have 



P(E™|E) 



^rri(Em)f^Af — m(E — Ej„) 



Taking the logarithm of both sides, we can expand 
f^w-m(E — Em) to first order in Em and find that 

InP(EmlE) = lnf7m(E„0 - V ^Ml^±^ (g) 
We define the canonical inverse temperature tensor by 



ainf7Ar(E) 



(9) 



and then the sum in equation [5] becomes the total con- 
traction aijE^ = Tr(Q;-^E„i)- The order of the indices 
is irrelevant since E is a symmetric tensor, and so is a 
through its definition. The canonical probability distri- 
bution for the total force moment tensor Em is then 

. . . -^Tr(<iS„) 

P^""(Em) = P(Em|E) - r!m(Em) (10) 



Z{a) 



The canonical partition function 



^m(a) = Y\ 

l,k>l 

acts as a normalization. 



(11) 



The angoricity defined by Edwards 4l| is related to 
the inverse of a. In its original definition, the angoricity 
is defined as the derivative of the stress tensor with re- 
spect to the entropy. From equation[21 a is the derivative 
of the entropy with respect to the force moment tensor. 
In the interest of controlling the proliferation of terms 
associated with temperature-like quantities for granular 
materials, we refer to the inverse of a as the angoricity. 
The modified Boltzmann factor for the granular system 
is then exp(— Tr(T~^Em))- Finally we show that the in- 
verse granular temperature in the canonical ensemble and 
in the microcanonical ensemble are equal. If we repeat 
the derivation of ([7]) using the form pO)) of the canonical 
distribution, then we obtain 



ainP(Ei|E) 



9 In r2 (El 



So the microcanonical inverse temperature equals the 
canonical inverse temperature, and is independent of the 
partitioning of the system. 

Special case of an isotropic system A lot of experi- 
mental and simulation effort has been devoted to systems 
under hydrostatic pressure [H, [13, El] caused by fixing 
the volume of the system in the absence of shear. Let 
r — Tr(E) be the trace of the force moment tensor. The 
extensive variable F is related to the intensive pressure 
by r = pA, where A is the area of the swtem. This 
makes F the internal virial of the system [43| . Then we 
can write the force moment tensor for an isotropic sys- 
tem in the absence of shear as E = ^/ [i^. It is natural 
to simplify the formalism to a single scalar variable F. 
Since the trace is additive, F is still a conserved, additive 
variable and the microcanonical and canonical ensemble 
derivations are the same as above. The key results are 
the definition of a 



51nriAr(F) 
dT 

and the form of the canonical distribution 



P^^"(Fm) = P(Fm|F) - r!m(F, 



Z{a) 



with the partition function 



Z{a) = J dFP'=™(F„ 



(12) 



(13) 



(14) 



We now show that a is related to the tensorial inverse 
temperature d by d = al: If we consider equation 
ini for an isotropic system, the density of states Sljv(E) 
must be symmetric under E^^ —>■ — E^^ since no di- 
rection of shear is preferable to another. Therefore, 
i7Ar(E) has an extremum at E^^ = 0, so that the log- 
arithmic derivative with respect to E^^ vanishes for the 
shear-free system and we obtain ai2 = (and by ex- 
tension a2i = since a is symmetric). Likewise, the 
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density of states must be invariant under rotations, so 
that the derivatives with respect to S^^ and are 
the same: an = 0122 — ol- Then the Bohzmann fac- 
tor exp(— Tr(dS)) becomes exp(— aF). The density of 
states for Fm, f2m(Fm), can be related to ^m^m) by 
using P^-"(F™) = Y.t^ P^^"(S™)5(F™ - Tr(i:™)) and 
we obtain 

f2m(r„0=n /rfS5^f^m(S™)<5(r™-Tr(S„)) (15) 

C. Discussion 

It is important to ask if there are important differences 
between angoricity and temperature. The most impor- 
tant distinction between a granular system and a thermal 
system is that the granular system has to be driven to 
change configurations. There is no simple equivalent to 
the thermal agitation which serves as a temperature bath 
for equilibrium systems, and which gives a natural value 
for the temperature. 

A granular system that is slowly sheared, so that 
it changes configurations based on the imposed strain, 
seems to come close to a thermal system. The exter- 
nal load resulting from the shearing sets the scale of the 
granular temperature, and the off-diagonal parts of 1/a 
can then be seen as a measure of the strength of the 
perturbations that the shear inflicts on the force chains. 
In this picture, the system stays in a force- and torque- 
balanced configuration until the load on a force chain 
becomes too large, upon which the system rearranges it- 
self into another configuration in mechanical equilibrium. 
Over time, the packing visits a large number of config- 
urations in a stress landscape analogous to the energy 
landscape for glassy systems. The dynamics of an equiv- 
alent system to the one described here, but in an energy 
landscape, is the subject of the SGR (Soft Glassy Rhe- 
ology) theory [3]. Recent work based on a toy model 
of activated dynamics in a stress landscape [i^ has been 
compared to experiments ^5^ as a test of the stress en- 
semble. The adaptation of the full SGR to the stress 
ensemble is in preparation [5l| . 

All of the above derivations can be performed in an 
equivalent manner for the volume as the extensive, con- 
served variable. The question of how the frequencies /3^^'^ 
with which the different microstates are accessed depends 
on the experimental or simulation protocol poses itself 
both for the stress ensemble and the Edwards ensemble. 

In the Edwards ensemble paradigm, this question was 
explored by O'Hern et al. [52| for very small disk pack- 
ings < 14 by enumerating all the states and mea- 
suring their frequencies for two different protocols. The 
same states were found with both methods, but as ex- 
pected with different frequencies. Thus the microcanon- 
ical equiprobability assumption is violated in this case, 
and there is no reason to expect a different result for the 
stress ensemble. 



We still have to investigate the validity of equation 
[5l the factorization of the density of states. This is a 
non-trivial assumption, since it breaks down if there are 
correlations between the subsystems we consider. For a 
volume-based ensemble derived along the same lines as 
above, Lechenault et al. have shown experimentally that 
even for sub-systems of size m > 100, there are corrections 
to equation [51 

A method to investigate the stress ensemble is sim- 
ilar to the approach taken for the Edwards ensemble. 
Where for the Edwards ensemble, the sytem is repeat- 
edly compactified at the same volume [10|, so that after 
each compression it enters a new configuration, we can 
create a new system with the same boundary stresses at 
each step. We now give a summary of our previous 
tests of the ensemble on simulated packings which are 
created from a random initial state and then relaxed un- 
til they reach mechanical equilibrium. 



D. Summary of simulation results 

We have tested the stress ensemble formalism on sim- 
ulated packings of frictionless disks with either harmonic 
or hertzian interactions in two dimensions using the 
algorithm of O'Hern et al. [l3|. We first rescaled the F- 
distributions of different configurations to test the form 
of the canonical distribution equation [T31 We find that 
we are able to perform the rescaling for any subsystem 
larger than to > 3, and that the equation of state a(F) we 
extract close to the unjamming transitiona (z) — Ziso = 4 
also becomes TO-independent above this value. From this 
we conclude that the factorization property equation [5] 
which underlies equation 1131 is valid for m > 3. 

For the systems with harmonic interactions, we have 
also fitted the density of states ri(F) to the form 

f](7)~F"° with a^2 + c{z- z,,af, (16) 

with c = 2.8 ± 0.5 (see Figure [1]). Together with the 
thermodynamic relation a — ^ 'gp this form leads to an 
equation of state over the full jammed range: 

«=^(2 + c((z)-z,,„)2) (17) 

The numerical density of states shows deviations from 
the form consistent with equation [T71 for m < 16. It is 
likely that this is a more sensitive test of equation [51 and 
the correlation length ^ « -\/l6 = 4 we can associate to 
this result is consistent with our results in the field theory 
section. 



E. Partition function at Ziso 

The results deduced from simulations, equations [TBI 
and I17[ depend on the specific force law, much like the 
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FIG. 1: Top: Distribution of F for a subsystem size m = 24 
out of TV = 4096, and fit using equation [16] for the density of 
states. Bottom: Fit coefficient a as a function of the contact 
number, with the fit which leads to equation 1161 



density of states of a solid depends on the detailed mi- 
croscopic interactions. Close to the jamming transition 
z — Ziso, however, we can derive the form of the density 
of states il(r) and the equation of state a(T) by count- 
ing the number of degrees of freedom only. This indicates 
that the coarse-grained properties of packings of friction- 
less spheres are universal, i.e independent of force law 
and simulation protocol, as Point J is approached. 

Counting the states If we use the definition of the 
density of states, equation 31 with equation [13] in the 
canonical partition function equation I14[ we are able to 
write 



Zia) 



/?^^"exp(-ar,), 



(18) 



where we sum over all packings consistent with force and 
torque balance, and with the force law respected (this 
last condition has to be modified for frictional packings) . 
We only consider packings of frictionless spheres, so the 
torque balance constraint is automatically satisfied. 



We can then formally separate the sum over config- 
urations v into a sum over all geometric configurations 
{rij} and a sum over all force configurations {Fij}, with 
(5-functions to enforce the force-balance (f.-b.) and the 
force-law (f.-l.) constraints: 



Z(a) 



EE' 



<M>'-'^-P-^,5(f.-b.)^(f.-l.) (19) 



For the frictionless packings, at the isostatic point, the 
number of degrees of freedom of the grains (dN) is equal 
to the number of forces that constrain them, N{z)/2. 
Equating these gives the isostatic contact number, Ziso — 
2d. It also shows that for a given geometry, there exists 
one and only one force configuration, and that the specific 
form of the force law becomes irrelevant. 

There is a one-to-one correspondence between a geo- 
metric configuration and a force configuration, so that 
we can eliminate the sum over the {r^ } and write: 



Z = 



J2 /3.'^"exp(-a^r,,({f^,,})F,,). 



(20) 



At the isostatic point, the mean force tends to zero, and 
the overlap (or deformation) of the grains becomes neg- 
ligible. Therefore, the dependence of {rij} on {Fij} can 
be neglected, ro, assuming that the system is 

monodisperse for simplicity. 

The simplest measure with which the space of 

states is sampled we can choose is of course the flat mea- 
sure /3jJ^" = 1. Although not necessarily correct, this 
simple approximation allows us to treat the problem an- 
alytically and we are able to extract the correct density 
of states (see below) . 

If we choose a flat measure, no interaction terms be- 
tween the different forces remain in equation 1201 Then 
the partition function can simply be written as a product: 



Nz 



oil 



z = 



n 

fc=i 



(iFfeCxp(-Q;roFfe) 



Nz 



o/2 



Connection to the simulation result The 



(21) 
partition 



function can be rewritten as a function of the scalar force 



moment F = pA'Y^f^^""^ riFi. We insert this equation 
as an identity into equation[2T]and obtain after switching 
the order of integration: 



Z{a) 



dFle""^ 
''o 



n 

i=l 



' NZi^a!^ 



dF,5 



E 



F 



We then perfom the force integral by enforcing the 5- 
function and then integrating over the other remaining 
forces with the constraint Fjq^.^^/2 G (0,F/ro): 



Z{a) = / d£- 

Jo ro 



dF. 
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The remaining integrations describe the volume of a piece 
of the Nziso/2 — 1-dimensional hypercube, with volume 

After dropping the prefactor the partition function is 
given by 

/■oo 

Z{a) = / r^-"o/2-ie-"r^ (22) 

This form is consistent with the simulation result for the 
density of states at the isostatic point (zi^o = 4 in 2d), 
n{T) = r^™, in the limit m >> 1. Either by using the 
thermodynamic relations a = ^ ^qy^^ on equation I22i or 
(r) = — ^glf on equation [211 one can also obtain the 
universal equation of state 



which matches equation [T71 for (z) Ziso- We have ob- 
tained the same equation of state for packings of disks 
with hertzian interactions close to Ziso |54| , showing that 
this result is independent of the interaction law. The de- 
viations from this form at small m that we observe the 
simulations confirm that the assumption (3^^^ = 1 breaks 
down at very small scales. 

However, the agreement we obtain at larger m shows 
that on a coarse-grained level, some of the properties of 
an isostatic packing can be understood through a sim- 
ple model which assumes a flat measure in configuration 
space. Counting the number of degrees of freedom is then 
sufficient to explain the form of the equation of state and 
the density of states (basically, the "thermodynamics" of 
the system). 

Distribution of the forces Equation [^T] also predicts 
the single-force distribution in the canonical ensemble 
if we assume a flat measure /JJ^"" — 1. From the 
form of the equation, we see that the probability to 
find a given force configuration is P{{Fi...Fj^^.^^/2}) = 
fl^""^^ exp(— aroFfc). Since this is a pure product dis- 
tribution, we deduce the single-force distribution 

P{F) ^ exp(-«roF) ^ exp ^- — ^j (24) 

i.e. a pure exponential. We have used the equation of 
state as well as the definition of the ensemble average of 
the forces {{F}) = {r)/Nro in the second equality. 

We emphasize that the canonical stress ensemble does 
not imply an exponential form for P{F) except at the iso- 
static point, if we assume a flat measure. An exponential 
distribution emerges at the isostatic point because for a 
flat measure the forces are independent random variables 
at this point. Similarily, in an ideal gas, the Boltzmann 
distribution, e~^^ , becomes an exponential distribution 
of energies because the energy E is a sum of single par- 
ticle energies. For interacting systems, e~^^ does not 



imply an exponential distribution of single particle ener- 
gies. 

Equation [24] is not a robust prediction since any devia- 
tion from the flat measure will have an especially strong 
effect on a single-particle quantity like the force distri- 
bution. With interactions, a calculation of P{F) is chal- 
lenging, as any one-body distribution is difficult to cal- 
culate for interacting systems [s^. The statistical me- 
chanical framework that we lay out in this paper is much 
more amenable to calculating correlation functions and 
response functions, and this is the task we focus on in 
this paper. 

A situation in granular materials that is very differ- 
ent from thermal systems is that numerically or exper- 
imentally, we do not have access to the canonical force 
distribution. Instead, by rescaling the force distributions 
by the spatial mean (F) = F/TVro of the forces (instead 
of the unknown ensemble average {{F))), we can mea- 
sure the microcanonical force distribution. O'Hern et 
al. established the algorithm used in our tests Q in p7j . 
where Figure 16 (top) shows the microcanonical force dis- 
tribution that can be obtained for the simulated system. 
It clearly decays faster than exponential. It is tempt- 
ing to use the results derived above for the flat measure 
and to translate them to the microcanonical ensemble. 
However, in doing so, we would neglect all correlations 
between individual forces which are clearly important at 
the single-grain level. 

Other theoretical approaches Kruyt and Rothenburg 
[56i] and Metzger et al. [H, H^l use a maximum-entropy 
approach with a multi-component Lagrange multiplier 
very similar to a to enforce that the total stress is con- 
served, and so work in the canonical a ensemble as well. 
The authors assume a product distribution for the forces 
and calculate the force distribution given the distribution 
of contact angles and distances between grains. The re- 
sult for the normal forces decays faster than exponential. 

Another approach to the problem is the force network 
ensemble (FNE) introduced by Snoeijer et al. [s^, which 
uses the decoupling of forces and the positions of the 
grains for very stiff grains. For a hyperstatic packing 
with (z) > Ziso, the FNE is then a microcanonical ensem- 
ble which assumes that for given a mean force (F) and 
a given geometry, all the conflgurations of positive com- 
patible with force and torque balance are equally likely. 

Tighe et al. [38!| simulate the FNE on the strongly 
hyperstatic (z) = 6 > Zi^o triangular lattice ad flnd a 
force distribution that decays faster than exponential for 
a system under isotropic compression, but an exponential 
decay for a sheared system. Recent work by Tighe et 
al. [53] introduces a second conserved quantity based on 
the height fleld to obtain a gaussian tail for an isotropic 
system. 

Snoeijer et al. [1^ derive an analytical force distri- 
bution in the FNE for an isotropically stressed trian- 
gular lattice, as well as for a general geometry. They 
obtain a density of states which scales as (F)^, where 
D N{{z) — Ziso) is the number of excess force degrees 



9 



of freedom in the system. 

Experimental results Measuring the force distribution 
inside a packing of grains is a challenge. Only two meth- 
ods have so far been successful: 

Majmudar et al. [l^, [2l| use quasi-twodimensional 
packings of photoelastic disks between cross polarizers 
and extract the stresses from optical measurements. For 
isotropic compression, they find that the distribution of 
the normal force components decays faster than expo- 
nentially, while the tangential force components follow 
an exponential distribution. If the system is subjected to 
pure shear, the distribution of the normal forces aquires 
an exponential tail while the tangential forces are not af- 
fected. The measurement is scaled by the spatial mean 
(Fn) of the normal forces, which relates this result to 
the microcanonical F-ensemble, as explained above in the 
context of the simulated data. 

Brujic et al. [60] as well as Zhou et al. [46] have mea- 
sured the interparticle forces using confocal microscopy 
on index matched suspensions of droplets. Again, the 
results are given scaled by the mean force in the config- 
uration. Brujic et al. have evidence for an exponential 
tail in the force distribution, while Zhou et al. focus on 
quantifying force chains. 

Experiments on quasistatically sheared systems in a 
Couette geometery have also produced force distributions 
[6ll . [13] , however the theoretical results above do not ap- 
ply to dynamical systems. 



III. BUILDING A FIELD THEORY 

By considering a field theory, we take a different route 
from the approach generally taken in the continuum me- 
chanics community. The focus there is to find a constitu- 
tive relation which links the stresses to the microscopic 
geometry of the packing. There is a considerable body of 
work on the subject in the mathematical and engineer- 
ing literature (see e.g. [1^ and references therein). The 
authors of [13] derive a constitutive relation for isostatic 
packings, though it can only be expressed at a microcopic 
level. 

A field theory, however, uses a path integral over 
all the possible configurations for stable packings. It 
coarse-grains the microscopic details of each packing into 
a continuous field which is sufficient to describe the 
macroscipic properties of that packing. Then we only 
need to combine symmetry arguments with a perturba- 
tive expansion in the flucutations of the field around its 
mean to obtain the weight of a configuration in the path 
integral. 

Here, we calculate correlations of the stress based on 
a minimal field theory that takes into account the essen- 
tial features of a frictionless granular packing. The field 
theory is dominated by a laplacian leading term which 
is multiplied by a stiffness which controls the behavior 
of the system as jamming is approached. We discuss the 
implications for the jamming transition. The field theory 



also predicts the correlation functions of the shear, which 
we test on simulated data. 

This field theory is related to our earlier proposal for a 
field theory [63] but, in this work, we have focussed only 
on stress correlations, and we have reassessed some of 
our assumptions based on information from simulations 
and experiments. The predictions from the earlier field 
theory related pressure to deviation from isostaticity, and 
its predictions for the jamming transistion have been fit 
to experimental data by Behringer et al. [2l|. 

A. The Airy stress function 

To write down a field theory of granular packings, we 
need to identify a field which incorporates as many of 
the constraints placed on the system as possible, such 
that they do not have to be imposed separately. The 
intuitive choice of the pressure p{x) — T{x)/A(x) as a 
field is misleading, since force and torque balance are not 
guaranteed for all possible configurations p{x). Instead, 
we use the Airy stress function 4', which incorporates 
force and torque balance constraints in two dimensions 
[12, [3, [13] , and which is related to the stress by 

a{x) = Vxh = VxVx*. (i.e.CTy = erkejidkdi'i'), (25) 

such that the pressure is given by p = Tr((T) = V^^*. The 
Airy stress function has been widely used in studying 
2D elasticity, and especially the role of defects [IB]- In 
its traditional usage, ^E" is obtained by minimizing the 
elastic energy. As will be seen from our analysis below, 
the field theory presented here has an entropic basis, and 
the role of ^ is very different. After expressing the path 
integral in function of the only remaining constraint is 
then that for purely repulsive granular packings, the local 
pressure has to be strictly positive. For three dimensional 
systems, has to be replaced by a 2-tensor *!' known as 
the Beltrami stress tensor Q|. We do not consider this 
case here. 



B. Minimal field theory 

We will work in the microcanonical ensemble, where 
the total stress F of the system and its contact number 
z are fixed. The key quantity to predict is therefore the 
microcanonical partition function Z{T, z) which is related 
to the canonical Z(a) by 

Z{a)^JdTdz Z{T,z)exp{-aT). (26) 

In a first step, we limit our investigations to two- 
dimensional isotropic packings under pure compression, 
such that the total force moment tensor S can be written 
as S = ^i. 

Let ip be the deviation of the Airy stress function 
from the one for a system with uniform pressure p = T/A. 
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Then the local stress tensor can be entirely written as a 
function of the second derivatives of ip 



—I+S(j = — 
2A 2A 



(27) 



From equation [21] we see that the Airy stress function 
admits a "gauge invariance" of the form 



■0(2;, y) ip(x, y)+ax + by + c; 



(28) 



that is we have the freedom to choose two arbitrary con- 
stants while constructing the microscopic of a packing: 
the position of the origin for the fluctuations of h and 
^p. This means that all physical quantities need to be 
at least second derivatives of 0. The field theory has 
to honor this symmetry and therefore can only contain 
terms with at least second order derivatives of ip. 

Here, we consider only systems under an isotropic com- 
pression, and hence the system has to be (statistically) 
isotropic. All the terms in the action have to honor this 
symmetry as well. Then we write (the T and z depen- 
dence is through the coefficients A and B) 



Z{T,z) = D^e^p 



dxdy ATi{Saf + BTT{Sa^) + . 

(29) 

In terms of the Airy stress function the 
two leading terms are given by A{W'^ijj)'^ + 

B ((9^0)2 + (9,^)2 + 2(a,a,v)')'. 

To extend this formalism to anisotropic systems, we 
need to decompose the stress tensor into a bulk term 
and the deviatoric stress Saf^^^ = Sdij — j^ijSakk- The 
coefficients multiplying the terms in the action are then 
analogous to the bulk and shear modulus of elasticity 
theory [65l |. 

The similarity of the above to the free energy formal- 
ism one writes in elasticity theory [3, EM] is due to the 
same combination of a tensorial quantity and symmetry 
arguments. The two are in fact fundamentally different: 
first of all, elasticity theory is written as a function of 
the strain from a given reference geometry. The elastic 
free energy is then used to determine the properties of the 
system if the reference geometry is disturbed. For a gran- 
ular system, there is no reference geometry and the strain 
is ill defined. Instead, equation sums over all possible 
geometries and forces at a given (Fjz). Moreover, while 
in an elastic system the Lame coefficients as well as the 
bulk and shear modulus are material constants, in this 
formalism they crucially depend on the imposed stress 
(see below). There is strong evidence that the bulk and 
shear modulus for granular and related systems depends 
on the imposed pressure and shear stress jl3| . 

In Fourier space, all the lowest-order terms are propor- 
tional to q^\ipq\^, so that we can condense A and B into 
one coefficient K and we write 



Z{T,z) = Jo^jexp 



^[K + A,q^+A,q^- 



(30) 



We have only investigated second-order correlation 
functions of the different components of the stress ten- 
sor in g-space. Since any eventual fourth order inter- 
action terms A(g)|'(/'g|^ (or higher) will just renormal- 
ize the Ak coefficients in the second-order correlation 
functions [53|, we are unable to probe them. Assuming 
possibly renormalized coefficients, we obtain a correla- 
tion function for the fluctuations of the local pressure 
6p = p- T/A = V^^: 



{\Sp,\ 



1 



Aq + A2q^ + AiQ'i 



(31) 



Simulation results 



We investigate the correlation functions S{q) of the 
local pressure pg = ^ J2ili ^gi^gi find out if they are 
of the type expected from equation [31] and to determine 
the dependence of the coefficients on F and z, as well as 
their interpretation. To obtain the data, we interpolate 
the discrete Pg onto a grid of size ^/N x ^/N and then take 
the two-dimensional FFT of the field. We then calculate 
the two-dimensional structure factor on the two- 

dimensional g-grid and finally take a radial average. 

The first observation is that the structure factor has 
an overall scaling form 



S{q) = T^s{q), 



(32) 



Figure[2]shows the limit K = limq^o S{q) from the lowest 
5 radial g-points (see figure [3]) as a function of F, with a 
F^^-scaling over 4 orders of magnitude. We only observe 
deviations in the limit of large F (and z) , far away from 
the jamming transition. To test if has the form 




FIG. 2: Parameter K for the g-independent term, determined 



from K = 1/ S{q). 
red. 



The scaling K ^ l/T is shown in 



equation [ST] we rescale by F^ and investigate the radial 
s{q), which is then only parametrized by z. We group 
configurations with similar z, and average over the s{q) 
to improve statistics (typically, about 20 configurations 
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are averaged over). With N = 1024, and a system size 
oi L X L grain diameters {L — 32 — 37, depending on 
packing fraction), we can investigate the range of q from 

2^ to '^^''"" x' A^grid = 32 is the size of the grid 

we use for Sp. 
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FIG. 3: Structure factor in log-log axes for systems with A'^ = 
1024 grains, grouped by z. The z increases from red to green, 
ranging from z = 4.02 for the upper curve to 2 = 5.72 for the 
lowest curve. 

Figure [3] shows the structure factor obtained for N = 
1024 for all z. The low-q values of s{q) decrease with in- 
creasing z, while the tail of the function does not change. 
The and g~'*-lines provided as a guide to the eye help 
to show the smooth transition of s{q) from a constant at 
low q through a q~^ decay at intermediate q to at 
high q. We were able to fit all the curves to equation!^ 
with 3 fitting parameters, and the results led us to write 
the following form: 



■<<!) = 



1 



fco(a(z) + e|g2 + ^4^4) 



(33) 
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FIG. 5: a) Fitting parameter K of the q-independent term, 
and fit to form equation 1331 b) Length scales ^2 and ^4 asso- 
ciated to the and g^-terms, in units of a. 



The following picture emerges: The two length scales 
^2 and ^4 are very nearly independent of the contact 
number z, and much smaller than the system size, with 
6 = 4.5 ± 0.5 and ^4 = 3.75 ± 0.25. This suggests that 
they reflect the purely microscopic properties of the sys- 
tem, like the distribution of grain sizes, which do not 
influence the behaviour of the system at larger scales. 
The q-independent term in the denominator, however, 
depends on z and is responsible for the supression of the 
\ow-q fluctuations with increasing contact number. The 
z-dependence is quadratic, 

k{z) ^ ko{2 + d{z - z,so)^) with 2 = 4±0.5, (34) 

with a coefficient of the same order of magnitude as the 
c of the equation of state [TT] 

The full second-order correlation function is then given 
by: 



Siq) 



1 



Ko (2 + 5(z-z,,o)2) + e|9'+e|9^ 



(35) 



FIG. 4; Fit of the structure factor to the form equation [33] 
for systems with A'^ — 1024 grains, grouped by z. 



D. Origin of the stiffness 



Figure [4] shows the fitted curves, while Figure [5] a and 
b show the parameters k = fcoa(z), ^2 and ^4. 



The previous section has shown that the behavior of 
the system at large length scales is dominated by the first 
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term in the field theoretical model, |V^VP- The length 
scales present in the system are just a few grain diam- 
eters, and remain small when the jamming transition is 
approached. There is no evidence of a growing static 
length scale in the system when the jamming transition 
is approached. This is consistent with our observations 
in the ensemble section, where we found that for length- 
scales larger than approximately 4 grain diameters, the 
density of states factorizes and the approximation of a 
flat measure in configuration space becomes valid. It is 
possible that there is a length scale unrelated to the pres- 
sure in the system, or a length scale that the finite size 
simulations are not sensitive to. Interestingly, the -0 field 
is critical and has power-law correlations, independent of 
the distance from jamming (see also next section). 

In the FNE, where by definition there are no correla- 
tions between the particles, it has been shown that the 
pressure correlations are flat, and that the correlations of 
the Airy stress function decay as I/9*, consistent with 
a field theory with only the stiffness K and no higher- 
order terms [66[. 

Our analysis raises many questions about the behav- 
ior of static granular packings near jamming (or more 
appropriately, unjamming). It is surprising that the sim- 
ulations show no crossover to usual elasticity theory, and 
that the 5* field remains critical even away from jamming. 
Usual elasticity theory would lead to a unique solution 
for ^E* at a given F and z. 

Can we understand why the form of the field theory 
is so very different from the common universality classes 
encountered in statistical mechanics? For example, in the 
(/)^-model the phase transition occurs when the mass term 
vanishes compared to the gradient term and higher-order 
terms [l5[. In our case, the mass term is absent because 
of the gauge invariance, and the tp-field is always critical. 

We need to understand the coefficient of (V^V)^: 



(36) 



and how it is related to the jamming transition at (F = 
0,z = 0). 

We can predict the scaling of the stiffness K{T, z) with 
F if we take into account the constraint that the local 
pressure p{r) has to be positive for all r. The argument 
below was originally proposed by B. Tighe in the context 
of the force network ensemble [66|. 

Let p(r) = T/A + V'^Sip > 0. After transforming into 
Fourier space, this condition becomes: 



(2^ 



< T/A 



(37) 



The LHS {— — V^V'(^)) is the negative of the local devia- 
tion from the mean pressure and can be both positive or 
negative. If the LHS is negative, the local pressure fluc- 
tuation is positive, and the constraint is automatically 
fulfilled. We therefore consider the case where the LHS 
is positive. We can square both sides while keeping the 



inequality, and transform one of the integrals by noting 
that Tp — -ip* 



F 



d^q 

(2^'' "^^^ i (2^'^ ""^^ - \A 

If we integrate over all r on both sides, the RHS aquires 
a volume term A, while in the LHS, we can change the 
order of integrations and get a (27r)^(5((f — q') from the 
exponentials. The condition becomes then: 

4| , |2 



q^\rP,\' < V'/A 



The integrand is always positive, so we can write: 



(38) 



The LHS is nothing but S{q) = l/K(r, z) and so, finally, 
the positivity constraint leads to 
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Kir,z)>- 



(39) 



The field theory we have constructed is for the marginal 
case where the stiffness satisfies the equality, and the stiff- 
ness is therefore the smallest allowed by the constraint 
of positivity. The z dependence of K is nontrivial and 
not predicted this argument. Its form, however, is not 
totally unexpected: At larger contact numbers, the stiff- 
ness, which is related to the number of configurations, 
should increase since there are more configurations avail- 
able to the packing. 



E. Tests of the field theory 

Implications for Jamming The form of the field the- 
ory we have obtained from simulations and from theoret- 
ical arguments presents a picture of the jamming transi- 
tion where the transition is the result of the number of 
possible states for the system tending to zero. As F tends 
to zero, the stiffness diverges, and hence the number of 
states around the smooth-pressure ground state that the 
system can access under perturbations is drastically re- 
duced. This is ultimately a consequence of the positivity 
constraint: the position of the hyperplane on which the 
F-constraint is satisfied shifts to the "lower left corner" 
of the space of allowed forces, so that the area of the 
hyperplane shrinks drastically. 

Fluctuations of the shear stress Even for isotropic 
systems, where the global shear stress Y^xy/A is zero, the 
local values of axy are nonzero and have well-defined cor- 
relations. Since the shear can be expressed as a function 
of the Airy stress function as axy = —dxdy'9, we can 
predict the shear structure factor in Fourier space: 



ili'^xy) 



2 9 

mi 



koq^ (2-Hc(z-z.,o)2) + e|<72 + ^|g4 ■ 



(40) 
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FIG. 6: Structure factor {\{(^xy)q'^) obtained from the mean 
of 20 simulated packings at (z) = 5.73. 




FIG. 7: Theoretical prediction of equation 1401 for {\{(Txy) 
at the same (z). 



Figure [6] shows the structure factor of the shear obtained 
from the mean of 20 simulated packings, while the the- 
oretical prediction from equation |40] is shown in Figure 

El 

We find good agreement between the simulation result 
and the prediction, especially for the angular structure 
of the correlation function. 

Real-space correlation Junctions From the form of the 
Fourier-space correlation functions equations [35] and BOl 
we can predict the real-space correlation functions e.g.. 



{5p{r)5pm 



(2^)2 



eyi^{iq.r){\5pq\ 



(41) 



and determine if the field admits long-range correlations. 
We condense the short correlation lengths ^2 and ^4 into 



a single correlation length ^ « 4, which multiplies the 
term. This is a good approximation since the integrals 
are dominated by the small-q-limit, where the q^-term 
dominates the q^-term in the denominator. 

Airy stress function The second-order correlation 
function of the fluctuations of the Airy stress can be cal- 
culated for ^ ^ 0, if the integral is cut off at system size 
(the cutoff is necessary even for ^ > 0). We obtain a 
scaling form in a; = ^j^. 



(^(r)V'(O)) 



ifo(2 + £(z-z,,o)2) 128(27r)3" 



[-^2^3(1, 1; 2, 3, 3; -x^) + 32 (ln(a;/2) + 2a;2+7-l)] 



(42) 



where 2^3(1, 1; 2, 3, 3; —x^) is a hypergeometric function 
and 7 is the Euler-Mascheroni constant. Figure [8] shows 
the correlation function for r G (0,i), and it is appar- 
ent that the correlations of the Airy stress function are 
long ranged and scale with the system size - qualitatively 
consistent with the real-space correlations seen in [sij . 

Pressure fluctuations The inverse Fourier transform 
equation BT] can be performed analytically. After the an- 
gular integration we obtain 



{5p{r)5pm 



dq 



Ka Jq 27r2 + c(z-z,^o)2+g2^2 



(43) 



The result is another Bessel function, Kq{x), the 0*'' 
modified Bessel function of the second kind: 



{Sp{r)6pm = 



Xo27r 



K.i^^±^^^^f^). (44) 



Asymptotically, for large arguments, ^^0(2:) ~ e~^ jx^^"^ 
so that the pressure fluctuations have short-range corre- 
lations which fall off beyond a scale set by the correlation 
length ^. 

Shear fluctuations For the shear fluctuations, the in- 
tegral has a nontrivial quadrupolar angular dependence 
and can only be done analytically for certain angular di- 
rections {6 = 0, 7r/2, TT and 37r/2) using polar coordinates. 
The result for these directions is a combination of Bessel 
functions and a constant piece, multiplied by a power 
law, 



{5d-xy{r)axyi0)) 



1 



+ 6 



f{z)r 



(45) 



where K^^ and are the second and third modified 
Bessel function of the second kind, and /(z) ^ 2 + 
c{z — ZisoY ■ The Bessel functions decay exponentially for 
r > ^, however the first term shows that the shear admits 
long range, power-law correlations ^ 1/?"^, regardless of 
the distance from the jamming transition. Equation 1451 
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has negative correlations, and in fact a numerical integra- 
tion of equation 201 shows that the real-space correlation 
function has a similar quadrupolar angular dependence 
as the Fourier-space form, with negative corrlations along 
the axes, and positive correlations at 45° (see Figure E] 
the oscillations in the correlation functions are the re- 
sult of a sharp upper cutoff of the g-integral at grain size 
a - ^) 




-4 -2 2 4 



X 



FIG. 8: Top: Real-space correlations of the Airy stress func- 
tion, expressed as a function of the scaling variable x = 
Bottom: Real-space correlations of the shear stress, obtained 
by numerical integration of equation 1401 



IV. MEAN FIELD THEORETICAL MODEL 

The results from the simulation, as well as the universal 
properties wc have found at the isostatic point can be 
combined into a mean- field theory of jammed packings of 
frictionless spheres. We investigate the properties of the 
minimum of the effective free energy F , and show that 
Point J has some properties of a critical point within this 
framework. 



A. An effective free energy 

At the isostatic point The results for the density of 
states from the simulation and the partition function at 
the isostatic point we derived in section 2 agree with 
each other (equations [IS] and . We can define an in- 
tensive mean- field variable 7 — T/m, such that ^(7) = 



^2mg amj _ |^^2g aj^"^ _ Then we Write a free energy as 
a function of this variable: 

Z{a) = / d7e-"-^('^) with 
Jo 

^^...0(7) =«7-21n(7). (46) 

Scaling of the mean contact number The variable 
which parametrizes the departure from the isostatic point 
is the mean contact number (z) = z, in mean- field no- 
tation. From the simulations, we were able to extract 
several scaling laws linking the contact number and 7, by 
exploring the phase space of compressed jammed configu- 
rations available to the conjuguate gradient minimization 
protocol. 

We observe the following relation between the ensem- 
ble means of (7) and (z) (see Figure IH]), 

(7) = B,{{z) - z^sof, with = 0.084 (47) 

{-i) ^ Bh{{z) - z,so? , with B,, = 0.026 (48) 

for harmonic and hcrtzian interactions, respectively. This 
scaling was first observed in the Chicago simulations [13] , 
from which our simulation protocol derives. 




FIG. 9: Coefficient for the power law scaling between (7) 
and (z) — Ziso for harmonic interactions (left) and hertzian 
interactions (right). 

This scaling, and by extension the empirical mean field 
free energy we present below, do depend significantly on 
the precise form of the intergranular potential (quadratic 
or hertzian in this case) and the nature of the simulation 
protocol. We explore the influence of the simulation pro- 
tocol in the context of a stability argument. 

Form of the free energy For z > Ziso, the number of 
additional force variables in a system of size N is N{z — 
Ziso)/2 and each geometric configuration can bear several 
force- balanced states. The form of the partition function 
at the isostatic point equation [21] needs to be updated 
to take account of the additional force variables, and the 
geometrical information cannot be dropped: 

Z = Y\. '^^^ exp(— 07) ^ 5 (geom.-constr.) (49) 
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The constraints imposed by the geometry on each ar- 
rangement of the forces are similar to the case studied 
by Snoeijer et al. [i^l- Their work focuses on integrating 
over all the possible force distributions allowed by force 
balance and purely repulsive forces on a given triangular 
lattice with fixed external forces. The case studied here 
is essentially an inversion of the problem: for a given set 
of forces, is there a geometry of {N, frictionless, spher- 
ical) particles to accommodate them? Their direct ap- 
proach, even if conceptually feasible for a single random 
geometrical configuration, fails for us because we have no 
Ansatz to tackle the problem of counting all the possible 
geometrical configurations at large compressions. 

We assume that the formal equation |49] is consistent 
with writing Z as a function of 7 and t — z — Zig^'- 

poo poo 

dj dt exp {mF^,^^ (7)) exp {mW{t, 7)) , 
Jo Jo 

(50) 

such that the effective free energy is F{j,t,a) = aj — 
2 In 7 -I- W{t, 7). We then determine the simplest W con- 
sistent with the simulation results discussed above. 

For a system with harmonic interactions, equation 1161 
directly leads to W{t) = — c(t)^ ln{'-f) + g{t), and this also 
gives the correct equation of state [T7] from setting the 
first derivative ^ to zero. To incorporate the relation 
BTl we set the ^-derivative of F to zero and substitute 
equation [47] for 7, so that we find 



g{t) = ct^ [HBst^) - 1] . 
Then, the effective free energy is given by 



_F(7, z) = a'f — 2 In 7 — ct^ 



In 



(51) 



(52) 



For systems with hertzian interactions, we only know 
the relation 1481 between 7 and z, but not the dependence 
of the density of states on z nor the equation of state. 
We can nevertheless obtain a similar free energy which 
incorporates equation 1481 



F(7, z) = a7 - 2 ln7 - Cht^ 



In 



7 



(53) 



This then makes a prediction for the density of states and 
the equation of state for a system with hertzian interac- 
tions: 



f](7, z) ~ 72+=-*' and a=—{2 + Ch{tf) (54) 

(7) 

More generally, for systems with a contact interaction 
of the type used in [l7| with a power S, we predict an 



effective free energy 

F(7,z) = a7-21n7-ct2(''-i) 



In 



/ 7 



VBi2(5-l) 



(55) 



B. Phase transition in the mean- field theory 

We now investigate the jamming transition in the con- 
text of the free energy equation [52l at first for a system 
with harmonic interactions. The constants in the free 
energy can be scaled out to give 



F(7,x) =a7~21n7-x2 In + 1 



(56) 



with the scaled variables x = c^^'^{z — ziso), 7 = -§7 and 



a = -a. 



Minimizing equation [56] with respect to its fields 7 and 
X allows us to extract the scaling of the ensemble averages 
(7) and (x) with the inverse temperature a. We find 



(7) = 



a — 1 



and 



a — 1 



1/2 



(57) 



and 



which is consistent with the equation of state at zi, 
the relation between 7 and z, as expected. 

Form of the transition The jamming transition is the 
point (7) ~* and (x) 0, i.e. a —^ 00. Figure [TU] shows 
the change in shape of the free energy approaching the 
jamming transition, and we observe a narrowing of the 
width of F around the minimum in the 7-direction, while 
the width in the x-direction increases. 
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FIG. 10: Contour plot of the free energy equation[56]for a = 5 
(left) and for a — 50 (right), approaching the transition. The 
scale for the contours in both cases is the height of the free 
energy barrier to a; = 0, F{{'y), (x)) — F{{'j),0). 

The signature of a second-order phase transition is a 
divergence of the fluctuations of the order parameters x 
or 7 as the transition is approached, i.e. a vanishing cur- 
vature of the free energy at its minimum at the transition 
[4^. We calculate the Hessian matrix at the minimum 
of F and we obtain that in the limit a —* 00, the eigen- 
vectors of this matrix become the 7- and x-directions, 
with eigenvalues that scale as ~ and = 4. This 
shows that the minimum in the 7-direction narrows dras- 
tically as the transition is approached, consistent with the 
well-definened equation of state we observe, and with the 
divergence of the stiffness in the field theoretical picture. 
The magnitude of the fluctuations in x, however, is in- 
dependent of the distance from the jamming transition, 
and does not diverge. 
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With 7 and x as order parameters, the jamming tran- 
sistion does not have the properties of a second-order 
phase transition. Neither can it be described as a first- 
order phase transition since it lacks a second minimum 
of the free energy. This reflects the fact that within the 
stress-ensemble framework, a granular packing can only 
be in a single state, i.e. jammed. Unstable configurations 
cannot be accounted for within the framework. 



stable 



• • • .*.• • *.*•• • *. • 




unstable 



C. A new order parameter 

Motivation The work of Wyart et al. [1, H^l and El- 
lenbroek et al. [4i] shows that there is an excess of low- 
energy deformation modes in soft granular packings close 
to the jamming transition. The authors link the soft 
modes to very low energy rearrangements in which the 
grains slip past each other. They also show that the scale 
of the rearrangements is determined by the distance of 
the packing from the isostatic point. A simplified version 
of the argument considers a packing of size I at mean con- 
tact number z, so that the number of force components in 
excess of the force and torque balance requirements in the 
packing is given hy Snj = NSz/2 ~ I'^Sz {5z = z—Ziso)- If 
we cut the boundary of the system, we remove of the or- 
der force components, and create an unstable region 
if Z"*"^ > I'^'Sz. This predicts a characteristic length scale 
I* ^ 5z~^ below which the packing is locally unstable, 
and this length scale diverges as the jamming transition 
is approached. The scale of the soft modes is set by T, 
so that their energy becomes uSl ^ (l*)^'^ ^ i^^)'^- 

The energy for the soft modes for a packing at pressure 
p = "f/A diminishes proportionally to the stress the soft 
modes cause at finite compression, which is proportional 
to p for a harmonic interaction potential. Therefore, the 
frequency of the soft modes at positive pressure scales as 
^ A{Sz)'^ — 7, where A is a constant. 

Every mechanically stable packing must satisfy the in- 
equality 



FIG. 11: Phase diagram derived from equation 1581 The red 
line marks the phase boundary between the stable and un- 
stable packings and the blue dots show stable packings. The 
sampling effect of the conjugate gradient minimization proto- 
col is illustrated by the black arrow leading from the initial 
configuration to the final stable packing on the phase bound- 
ary. 



the contact number z will increase. Stopping the pro- 
cess at the first stable configuration reached biases the 
process towards configurations with high 7 and low z 
compared to a flat sample of all the stable conflgura- 
tions. The same bias is not immediately encountered 
in other methods used in the literature to create me- 
chanically stable conflgurations, such as tapping [l3| or 
shearing [so'l . More generally, the marginally stable con- 
figurations should be encountered in any protocol which 
does not allow the system to thermalize, i.e. explore the 
phase space of jammed configurations, but instead does 
an infinitely rapid quench to the nearest stable packing. 

Figure [H] shows that as we approach the isostatic point, 
the relative fluctuations around the stability line derived 

above increase. Therefore, it is natural to investigate the 

2 

variable u — which takes the value of 1 on the stability 
line, and u > 1 in the stable region. 

Mean field theory in u-x coordinates We can rewrite 
the free energy equation [51] in terms of the new variables 
u and X 



0<6E < A{z^ 



1, 



(58) 



where 5E is the energy of the lowest-energy displacement 
eigenmode of the system. 

The equation predicts a phase diagram for static gran- 
ular packings, with stable packings above the 7 — A{z — 
^iso)^-line, and unstable packings below it (see Figure 

[ni. 

The simulated packings always lie on the boundary 
between the stable (7 < B{z — ZisoY) and unstable 
(7 > B{z — ZisoY' regimes, i.e they are marginally sta- 
ble. Wyart et al. have suggested a mechanism to ex- 
plain this 0: The conjuguate gradient minimization will 
decrease the total potential energy of the system by re- 
ducing the overlap between initially random disks just 
enough to reach mechanical equilibrium. During the pro- 
cedure, on average, 7 will reduce with the energy, while 



F{u,x) 



2 In 



(1 — Inu) 



(59) 



and investigate the properties of the transition in func- 
tion of u and x. The position of the minimum is at 

(m) — 1 and (x) — I J , which is compatible with 

the results obtained in the (7, a;)-coordinates. Figure fT2l 
shows the evolution of the shape of F as the jamming 
transition is approached. We observe a striking asym- 
metry develop in the width of the minimum of F as the 
jamming transition is approached. It is clear that there 
exists a direction along which the second derivative of F 
vanishes, and hence the fluctuations diverge. The hessian 
matrix is given by 



d^F / 2^ 
du,dx l-2[2(a-l)] 



1/2 



-2[2(a- 1)]'/'^ 
4(a~l) , 



(60) 
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FIG. 12: Contour plot of the free energy equation[59]for a = 5 
(left) and for a = 50 (right), approaching the transition. The 
scale for the contours in both cases is the value of the second 
derivative of F with respect to u at the minimum of F. 

The eigenvalues close to unjamming, for large a, be- 
come A- ~ ^ and A+ ~ 2a. Hence the curvature of 
the free energy vanishes along the direction of V- , whose 
angle with the M-axis is just tan(6') — v'^/v^. At the 
unjamming transition a —^ oo, this angle vanishes as 
ta,n{9) ~ {2a)^^/^, and the flat eigendirection becomes 
the u-axis. 

So in the limit a — > oo, we observe a diverging suscep- 
tibility in the u-dircction, that is the fluctuations around 
the stability line 'y — x'^ diverge. This indicates criti- 
cal behavior at the jamming transition, further evidence 
that Point J is a critical point if approached along the 
T = 0-hnc. 

Comparison to simulation data The first prediction 
of the field theory in u — x-coordinates is that all the 
(m, x) data points associated to individual configurations 
should cluster around u = 1, after an appropriate rescal- 
ing. Second, we predict that the x value of of the data 

/ \ 1/2 

points for a given a clusters around {x) — ( j . Fig- 

ure[T3]shows configurations from throughout the jammed 
region, grouped by their values of a (indicated by the 
color of the data points). The data cluster around u = I, 
with large fluctuations of for data points at higer a (in 
red). The circles, which are in the same color as the data 
points mark the minimum of F for the a associated to 
that color. The line through each circle is in the direction 
of the eigenvector along which the susceptibility diverges. 
Its length is proportional to 1/A_, i.e. proportional to 
the magnitude of the expected fluctuations. We observe 
that the data points do indeed cluster around the min- 
imum and follow the direction given by the eigenvector. 
In the limit of the jamming transition, the spread of the 
data points at a given u becomes very large, and parallel 
to the M-dircction, as expected. 

We can obtain a scaling form of the free energy close 
to the jamming transition. The a;^ (1 — In w) term in equa- 
tion [55] becomes vanishingly small compared to the other 
terms in this limit, and if we define t = a^l'^x we can 
write the free energy (up to a constant) as 

t^ 

i?(i,u) = _-21n - . (61) 
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FIG. 13: Data points from the A'^ = 1024-configurations with 
harmonic interactions grouped by different values of a (red: 
high a, green: low a) plotted in the x — u diagram. For each 
a, the circle of the same color marks the minimum of F at 
that a. The line through each circle is in the direction of the 
eigenvector along which the susceptibility diverges. Its length 
is proportional to 1/A_, i.e. proportional to the expeted fluc- 
tuations. 
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FIG. 14: Same as Figure [13] after the a;-axis has been rescaled 
by a^^^ , superimposed on a contour plot of the scaling limit 
of the free energy equation |611 (see text). 



Figure [14] shows a contour plot of F{t,u), and superim- 
posed on it the same data points as in Figure [T31 mul- 
tiplied by their individual a^^^. We obtain an excellent 
collapse of the data for all the configurations near Point 
J, and the divergent fluctuations along the minimum of 
F can clearly be seen. 
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D. Link of the mean field theory to the field theory 

For the field theory derived in section 3 and the mean 
field theory explored above to be consistent with each 
oher, the partition function has to take the following 
form: 



Zm{a) — / d'-fdzexp{—mF{j,x)) 



d'^dz / Dip exp 



2 



+^2V+e4V)i^,i 



r 

exp( 



2tt' 



((2 + c(z - z,sof)q^ 
(62) 



We can show this link in two ways. The simpler method 
is to expand the mean field theory to second order about 
its minimum in 7(7^ = 7 + ^'^ip{r), while keeping the 
scaled contact number x as a parameter. We obtain 



7 



(63) 



which is nothing but the stifi^ness part of the field the- 
ory. This shows consistency since the purely microscopic 
length scales ^2 and ^4 are not expected to contribute to 
the mean field free energy. 

Alternatively, we can perform the path integral in ip 
and see if we obtain the entropic part of the mean field 
free energy from the logarithm of the microcanonical par- 
tition function. Although this can be done exactly, since 
we limited ourselves to gaussian terms, there are prob- 
lems with this approach. The field theory is an expan- 
sion in powers of ip, and we have kept only the lowest 
order. This is sufficient to calculate the moments of ip, 
like the structure factor, but it is probably inadequate to 
accurately calculate the generating functional Z{'-f, x). A 
proper treatment of the problem would need knowledge 
of the full power series and then a sophisticated renor- 
malization approach, for a two-dimensional problem. 

Nevertheless, some progress can be made and we can 
test the consistency between the two forms to lowest or- 
der in 7 and we obtain [53| 

^(7, z) = ^0 + (2 In 7 - ln(2 + £(z - z,,of)) . (64) 

Comparing this to the mean-field result, we find that 
the isostatic limit of both expressions, 2 In 7, is identical. 
We recover part of the ^-dependence if we expand the 
logarithm of the contact number ln(2 -t- c{z — Ziso)'^) ~ 
In 2 -I- 1/25(2; — Ziso)'^- The field theory does surprisingly 
well in a numerical comparison to the mean-field result. 
Figure [T5] shows both for the pairs of (7, z) corresponding 
to the structure factor plotted in Figure [31 This is likely 
due to the fact that the (7, 2:)-pairs investigated are all 
quite close to the line of u = 1, so that the logarithmic 



term In 



vanishes. On the line u = 1, both 



expressions reduce at first order to 

S ^ So + N \2lii-/ - cx^ 



with a value c = 2.8 ±0.5 from the mean field theory and 
a value c — c/2 — 2± 0.25 from the field theory. 




(65) 



FIG. 15: Comparison between the mean-field theory (blue 
line) and the field theoretical result for S(7, 2), using the val- 
ues of (7, z) from the simulations for the latter (red dots). 



V. CONCLUSIONS 

A. Summary 

We have investigated the properties of jammed gran- 
ular assemblies, approaching Point J and within the 
jammed region, using a newly formulated stress ensem- 
ble. The jamming transition can be analyzed within this 
framework. 

The force moment tensor of a system can be written 
as a boundary term, which makes it a conserved quantity 
under internal rearrangements. We use this conservation 
law to define a canonical ensemble in section [Til The 
conjuguate variable to the force moment tensor defines 
then a granular analog to temperature, that is closely 
related to the Edwards definition of angoricity. 

We review our previous test 0] of the ensemble, and 
obtain an equation of state for the granular temperature 
within the jammed region. The form of the equation 
of state is a universal property at the isostatic point, a 
conjecture that we confirm through an exact calculation, 
though we also detect evidence for short-range correla- 
tions. 

In section 3, using the constraint-free Airy stress func- 
tion, we build a field theoretical model for the jammed 
range based on symmetry arguments and the analysis 
of the pressure fluctuations in the simulated configura- 
tions. The jamming transition appears as a divergence 
of the stiffness of the leading term, so that the entropy 
of jammed packings tends to zero at the transition. 

Finally, in section 4 we combine all the relations ob- 
tained from simulation into a phenomenological mean- 
field theory. We investigate the jamming transition again 
to determine the correct order parameter. Inspired by a 
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stability argument proposed by [1], we use an order pa- 
rameter which measures the deviatins from the stabihty 
Une Unking mean pressure and mean number of contacts. 
The mean field theory then predicts divergent fluctua- 
tions in this order parameter at Point J, lending weight 
to the interpretation of J as a critical point. 



B. Limitations of the ensemble, predictions and 
further tests 

Limitations and applicability The stress ensemble, 
with its tensorial conserved quantity and tensorial tem- 
perature, can be applied to a wide range of systems. The 
only conditions are that the system has to be in mechan- 
ical equilibrium, and that there can be no long-range in- 
teractions between the particles. This excludes any kind 
of driven system (exceptions below), or system with a 
temperature, but it includes systems with attractions, 
like colloids, or packings with friction. Anisotropic sys- 
tems, either made of nonspherical particles and/or with 
an imposed shear in addition to a hydrostatic pressure, 
are a major area where we expect the stress ensemble to 
be useful. 

The analysis of the simulated packings, the field the- 
ory, the mean field theory and the interpretation of Point 
J in this context are limited to isotropic, frictionless pack- 
ings of round grains. Lifting any one of these restrictions 
fundamentally changes the nature of the system. 

For packings with an imposed shear, we expect the 
"Boltzmann factor" exp(— Tr(Q;E)) to reduce to a factor 
exp(— ttpp) featuring the presssurep = 1/2(cti -1-0-2) and a 
pressure-temperature ap and a factor exp(— agr) featur- 
ing the shear t = 1/2 (di — (T2) and a shear-temperature 
as- Here (Ti and (T2 are the principal stresses for the 
global stress tensor. We have extended the field theory 
presented here to the pure shear case, and are in the 
process of testing its predictions against simulations and 
experiments (65j . 

The introduction of friction, or of anisotropic grains 
transform the isostatic point into a broader region within 
which marginally stable packings occur. It is then doubt- 
ful if Point J is a single, well defined singular point for 
these systems, and it will be interesting to investigate the 
field theoretic framework for these cases. 



C. Predictions 

In principle, all the relations obtained from the ap- 
plication of the ensemble to the simulated system are 
predictions which can be verified in other simulations or 
in experimental systems which are close to frictionless, 
round particles, like bubble rafts. There are however sev- 
eral major caveats. 

First, the density of states at larger compressions is 
dependent on the method of preparation of the sample, 
as discussed in section 4 and it remains to be seen which 



range of (simulation or experimental) protocols selects 
for marginally stable packings. The behavior close to the 
jamming transition should not be affected, though. 

Second, most of the predictions were made in the a- 
canonical ensemble which is difficult to reproduce in ex- 
periments or in simulations since there are no means of 
imposing a externally, as is normally done with the tem- 
perature in thermal systems. It is feasible, though, to 
work in the microcanonical ensemble where the external 
stress (or the hydrostatic pressure, for isotropic systems) 
is fixed. The predictions of the canonical ensemble can 
be adapted to the microcanonical ensemble by carrying 
out a Legendre transform if the subsystem investigated 
is sufficiently large (i.e. m » 1). For smaller subsys- 
tems, especially for the force distribution and the single- 
grain pressure distribution, the local correlations make 
the thermodynamic approximation break down. 

The best prediction is the equation of state a = 
^'°P^ , together with the coarse-grained density of states 
^iXm) ~ r^™ close to Point J for 2d frictionless round 
grains. It is a universal property which we expect to hold 
regardless of preparation method and the choice of the 
ensemble. 



D. Further tests 

There are several interesting methods by which the 
stress ensemble can be explored. 

When a hot and a cold body come into touch, the sec- 
ond law of thermodynamics dictates that heat will flow 
from the hot body to the cold body until both of them 
arrive at the same temperature, regardless of the com- 
position of the two bodies. In the context of the stress 
ensemble this means that two compartments at different 
cki ^ a2, if brought into touch, will transmit stress from 
the compartment with lower a (i.e. higher granular tem- 
perature) to the compartment with the higher a, until 
Qfi = a2- To set up a numerical experiment testing this 
situation, one could prepare two sets of simulated con- 
figurations similar to the ones tested in determine 
their respective equations of state, and then bring them 
into touch. Then one can measure the Fm-distributions 
in the two halves, determine their respective a via the 
equations of state and test if both values of a match. 

The only dynamical situation where the stress ensem- 
ble can be applied is a quasistatic motion where the sys- 
tem evolves through a sequence of equilibrium states. 
One example of this is a system under quasistatic shear, 
slow enough for it to be characterized by a sequence of 
stress buildups and eventual rearrangements. It is in this 
regime that we expect that the SGR formalism [l3| can 
be adaped from the Boltzmann ensemble to the stress 
ensemble, as we have done for a toy modelji^. 

Another situation of the same type is a system which 
is periodically tapped so that it rearranges itself into a 
new configuration. In contrast to the same setup in the 
context of the Edwards ensemble, we need a system with 
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the same externally imposed stress after each rearrange- 
ment. For this system, or for the quasistatically sheared 
system, we can then employ a granular version of the 
fluctuation-dissipation theorem (FDT), as has been done 
by Song et al. for the Edwards ensemble [35i|. Let F be 
an external perturbing force, and x{t) be the position of 
a tracer particle. Then the FDT links the fluctuations in 
the particle position to its mean response to the pertur- 
bation: 

(WO-(0)f> = -^^^^^^^ (66) 
a t 

The driving force, provided by gravity in the work of Song 
et al., should be replaced by e.g. the magnetic force on a 
metallic tracer particle since the stress ensemble is sen- 
sitive to a gravitational field. The values of a extracted 
by this method should then only depend on the exter- 
nal stress, not on the type of material employed or the 
magnitude of the driving force. 

E. Outlook 

This work underscores that the experimental condi- 
tions under which a granular material is examined are 



crucially important. The framework we have introduced 
permits us to distinguish between systems in a canonical 
stress ensemble or in a microcanonical stress ensemble, 
and between different stress states. The statistical frame- 
work can be used to predict correlations functions of the 
stress under these different conditions. 

The jamming phase diagram masks a more compli- 
cated reality, arising from the sensitivity of the jamming 
transition to the prepared state of the packing (see @ for 
a review) . The stress ensemble provides a concrete frame- 
work for understanding the behavior of granular materi- 
als at the jamming transition and within the jammed 
region. 
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